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Abstract. A fundamental property of any material is its response to a lo- 
calized stress applied at a boundary. For granular materials consisting of 
hard, cohesionless particles, not even the general form of the stress response is 
known. Directed force chain networks (DFCNs) provide a theoretical frame- 
work for addressing this issue, and analysis of simplified DFCN models reveal 
both rich mathematical structure and surprising properties. We review some 
basic elements of DFCN models and present a class of homogeneous solutions 
for cases in which force chains are restricted to lie on a discrete set of directions. 

1. Introduction. Cohesionless granular materials exhibit a range of behavior that 
has fascinated physicists, applied mathematicians, and engineers from Coulomb, in 
the late IS"* century, to Schaeffer, in the late 20*'' and early 21**. It is therefore 
somewhat surprising that perhaps the simplest question one can ask about a gran- 
ular material remains unanswered even now. Given a box full of sand, if a marble 
is placed on top of the sand, how is the weight of the marble distributed on the 
bottom of the box? The problem is not just quantitative; we do not even know 
the qualitative form of the pattern of vertical force on the bottom generated by the 
marble. We do not know, for example, whether the vertical force is largest directly 
underneath the marble or has a maximum on a ring. ^ 

Experiments suggest that qualitative features of the response may depend on 
the way in which the sand was put into the box or on the geometric and frictional 
properties of the grains. [TH Experiments on 2D layers with depths up to 

approximately 12 grain diameters show both two-peak and single-peak responses, 
depending on the geometry of the grains and/or degree of disorder in the packing. 

Experiments on 3D layers have probed the response in layers as deep as 300 
grain diameters, and single peaks directly beneath the applied load were observed. 
[TH] A theoretical framework that allows a unified interpretation of these results 
would be of great interest. 

Let us formulate a thought experiment that focuses on the central features gov- 
erning the stress pattern at the bottom of the sandbox. The fundamental problem 
has to do with how stress is transmitted through the bulk of the material, so it is 
natural to eliminate the side walls from consideration by imagining a slab of sand 
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^ A review of the many and varied approaches that have been taken to the problem of determin- 
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that is infinite in horizontal extent. It is also tempting to eliminate the effect of 
the bottom boundary by imagining it to be at infinity, where all stresses are in- 
finitesimal, and calculating the stress as a function of depth far from the bottom 
boundary. We will see, however, that in the context of an interesting class of models 
the bottom boundary condition plays an essential role. 

We can also eliminate from the problem the complicating effects of gravity. We 
consider a slab of sand between two plates that are being pushed together. The 
density of the sand is assumed to be high enough that the material rigidly resists 
the force applied by the plates. The material is then subjected to a localized force 
created by poking a needle through the top plate and applying a specified force 
to it. Note that both plates and the confining pressure are essential in this setup. 
Without them, some grains at the surfaces of the slab would likely be subject to 
outward forces that could not be balanced, so the material would not support any 
force at all on the needle. 

A third simplification is to assume that the bulk granular material is isotropic. 
That is, we assume that the geometry of the grain packing does not distinguish any 
particular direction in space. While the local environment of any particular grain is 
clearly anisotropic, we assume that this anisotropy, as measured for example by the 
fabric tensor which characterizes the degree to which there are preferred directions 
for "bonds" joining the centers of grains in contact, vanishes when coarse-grained 
over sufficiently large volumes. This is a strong assumption, as construction history 
of a real granular material may well distinguish favored directions for contacts. |16| 
Nevertheless, there are highly nontrivial structures to be uncovered even in isotropic 
models and it is surely appropriate to understand them before attempting to include 
the complications and additional parameters associated with intrinsic anisotropics. 

Finally, for simplicity, we can work with a two-dimensional system. The top of 
the slab is a line defined to be at z = and the bottom is defined to be z = d. The 
slab is infinite in the horizontal x direction. We are interested in determining the 
response to a localized force applied at x = and z = 0. 

Treating the sand as an ordinary (linear) elastic medium, a standard treatment 
indicates that azz{x,z) should have a single peak centered on a; = with a 
width that grows linearly with depth. ^ This treatment assumes, however, that the 
stress field is governed by equations derived using the continuity of a well defined 
displacement field and that there is a unique stress associated with each strain 
configuration. Neither of these is true for a granular material consisting of rigid 
particles that support frictional forces, so a new approach is needed. 

The problem of deriving macroscopic stress equations from known microscopic 
or grain scale physics has proven quite difficult. An indication of just how difficult 
this might be can be found in almost any experiment or numerical simulation that 
generates images of the intensity distribution of stresses on scales of several grain 
diameters. ^1 ^| In such images, one sees immediately that the stress is 
concentrated on filamentary structures called force chains that support compres- 
sive stress. These chains appear to be relatively straight on scales of up to 10 grain 
diameters or so and give the visual impression of splitting and fusing at a variety 
of angles. The presence of such structures suggests that passage from the grain 
scale to the macroscopic stress will involve two distinct steps: we must understand 



^In sufficiently strongly anisotropic elastic materials, it is possible to have two peaks, both of 
which have widths that grow linearly with depth. 1171 
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how grain scale physics favors the formation of chains, and also how the interac- 
tions among chains determine the macroscopic stress field. While these two tasks 
must ultimately be facets of a unified theory, it may be useful to approach them 
separately. 

Models of directed force chain networks (DFCNs) have recently been proposed 
to bridge the gap between the scale of individual chains and the macroscopic stress, 
leaving open the issue of how grain scale physics promotes chain formation. |2UI l?T| 
(This type of model is also referred to as a YA-model.|2n|)- As detailed below, 
a "Boltzmann equation" governing the densities of chains with specified strengths 
and orientations can be obtained, assuming only that whether a chain splits at a 
given point is determined only by the local environment around that point and that 
environments that lead to splitting are homogeneously distributed throughout the 
material. An essential feature of this equation is that the unstressed solution is 
unstable; perturbations of it lead to divergent responses. Thus in order to calculate 
physically meaningful response functions, it is necessary to perturb around a pre- 
stressed state. The only case for which response functions have been calculated is 
the special one in which all force chains support the same magnitude of stress and lie 
on a 6- fold symmetric set of vectors. In that case, the response consists of two peaks 
whose centers diverge linearly with depth but whose widths grow like the square 
root of depth. In other words, the only solved case of response functions for DFCNs 
suggests that stress propagates along characteristic directions determined by the 
homogeneously pre-stressed state, in marked contrast to expectations from standard 
elasticity theory. This propagation of stress is a feature of noisy hyperbolic 
equations and is therefore called a "hyperbolic response." [22] 

It is worth noting here that hyperbolic response is known to occur in models 
of isostatic packings of frictionless grains and that simple constitutive rela- 
tions leading to hyperbolic stress equilibrium equations have been shown to explain 
nontrivial features of experiments on sandpiles and granular columns. |23| One 
must also note, however, that direct measurements of the response seem to indicate 
that single-peaked response is associated with strongly disordered packings, and 
hyperbolic response occurs only in ordered systems or at shallow depths. |14l 

The question immediately arises as to whether the hyperbolic response of the 
6-fold DFCN model is an artifact of the confinement to a discrete set of directions, 
or perhaps of the unique property that all chains in the 6-fold DFCN have the same 
strength. To investigate the robustness of the 6-fold DFCN results, one would like 
to find other analytically tractable cases. Here we address the first step in this 
direction - the identification of models for which stable solutions can be found for 
homogeneous macroscopic stresses. 

In this paper we present homogeneous solutions to the DFCN Boltzmann equa- 
tion for special cases in which chains lie along a discrete set of directions with 8-fold, 
or more generally 4N-fold, symmetry in two dimensions. Even the homogeneous 
solutions - the pre-stressed states about which we might contemplate calculating 
response functions - have nontrivial features. The actual calculation of the response 
functions these solutions is beyond our present scope. 

This paper also provides corrections to the presentation fo the general Boltzmann 
equation of Ref. [23 • The corrections do not affect any of the calculations on 
discrete models reported in that paper, but do clarify some conceptual points crucial 
for future work on continuum models or models permitting multiple splitting angles. 
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Figure 1. Illustration of a directed force chain network. Line 
thicknesses indicate force chain intensities. Arrows indicate force 
chain directions. Large circles mark regions of the granular mate- 
rial containing local configurations that require the splitting of a 
force chain entering from the direction pictured. Small circles in- 
dicate local configurations where fusion of two incoming chains in 
the directions pictured occurs. The boundary forces applied to the 
system are shown as solid vectors above the top surface and below 
the bottom surface. The dashed vectors at the two surfaces indi- 
cate forces exerted by the top and bottom plates on the material 
inside as a response to the applied forces. 

2. The DFCN model and Boltzmann equation. 

2.1. Definitions. A DFCN is a collection of line segments with associated com- 
pressive forces assigned such that force balance is achieved at every vertex. Figure^] 
shows a simple example that illustrates several features. First, each segment is as- 
signed a strength, or force intensity, that is indicated in the figure by the thickness 
of the line. Every chain is assumed to be under compression, so the force exerted 
on a vertex due to a given chain is directed along the chain toward the vertex, and 
the vector sum of the forces at any given vertex must vanish. 

Second, each segment is assigned a direction, indicated by the arrow, which 
distinguishes its "beginning" from its "end." The physical distinction lies in the 
local grain configurations at either end of the chain. Certain configurations do not 
allow propagation of a chain through them in certain orientations and therefore give 
rise to force chain splitting. The relevant configurations of this type for the network 
shown in the figure are indicated with large circles. Other local configurations would 
not necessarily require a splitting event, but do permit two chains to fuse when they 
intersect. The relevant configurations of this type are indicated with small circles. 

The directions of chains in the entire network are determined by the specification 
of the applied boundary forces. One must think, for now, of the force being applied 
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Figure 2. Illustration of force chains, (a) Two force chains ini- 
tiated at the top surface cross without interacting, (b) When one 
of the force chains is initiated at the bottom, a fusion becomes 
possible, (c) A force chain initiated at the top that splits upon 
encountering a defect. 



at point A as coming from a needle being poked through the top plate and held 
at constant force. This ensures that a chain beginning at A must be present. At 
point B, a local configuration exists that requires the chain to end. To balance the 
force at B, two chains must be initiated. Hence the arrows on these to chains must 
point away from B. 

By tracing the chains initiated by applied forces at the boundary and taking into 
account the type of local environments at the chain endpoints, each chain direction 
can be uniquely assigned, with the exception of rare cases in which fusions happen to 
occur precisely at places that might be mistaken for splittings. When all directions 
are assigned, one is likely to find some chains that end, rather than begin, on a 
boundary. The boundary forces required to balance these forces, shown as dashed 
arrows in Figure must be interpreted as a response to the applied forces. 

Finally, there is the possibility that a single chain can appear to have inconsistent 
requirements on its direction, as occurs for the chain passing through point E in 
Figure n Here two splitting events have initiated chains along the same direction 
that meet head on. Formally, this corresponds to a fusion event in which the 
outgoing chain has zero strength, which thus appears as the annihilation of two 
chains, and the circle marking this fusion can be imagined to lie anywhere along 
the chain. The probability of such configurations is determined by the grain size, 
or, more precisely, by the ratio of the grain size to the typical separation between 
force chains of opposite directions. In the following, we will take this probability 
to be zero. (But keep in mind that a nonzero probability can be invoked to cut off 
the divergence in the density of weak force chains that arises in the 8- fold model.) 

FigureElfurther illustrates the difference between chain splitting and chain fusion 
in the context of a hexagonal array of grains. In panel (a), two forces are applied 
at the top boundary and the forces propagate along chains, crossing at the central 
grain. In panel (b), a similar situation is shown, but this time one of the chains is 
assumed to be specified by fixing its position at the bottom boundary instead of the 
top. In this fusion of the two chains at the central grain is permitted (though 

not required) and the resulting configuration may be different. Panel (c) shows a 
local configuration (a missing disk) that requires the splitting of an incoming chain. 
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Thus each chain is characterized by an intensity, /, and a direction n, or, equiv- 
alently, a vector /. Force balance at a vertex requires that ^ / for incoming chains 
equals ^ / for outgoing chains. Positive values of / correspond to compressive 
stress along a chain, negative values to tensile stress. To model noncohesive mate- 
rials, which do not support tensile stress, we take / to be always positive definite. 
Thus there is no ambiguity introduced by using / to denote the pair (/,«). When 
convenient, we will use the angular variable 9 to indicate the direction of fi. We will 
use the term "/-chain" to refer to a chain with the given strength and direction. 

2.2. The Boltzmann equation. The following discussion supersedes the treat- 
ment given in Ref. |21j. In that paper, mistaken reasoning was used in writing the 
continuum equation. The mistake does not affect the bulk of that paper, since the 
forms obtained there for networks of chains confined to a discrete set of directions 
were the same as those derived below. Nevertheless, the corrections made here may 
be important for future work on continuum models and are necessary for conceptual 
clarity. Note especially that the definition of P has different dimensions here from 
that used in Ref. [21] • The theory is developed here explicitly for the case of two 
dimensions. Generalization to higher dimensions is straightforward. 

Let P{f, 9, r) represent the density (in an ensemble average) of force chains of 
intensity / and direction 9 passing through the spatial point r. In other words, 

/ f'df d9' \n-u\P{f\9' (1) 
Jf Je 

is the number of chains with intensity between / and f + Sf that cross a unit length 
line segment passing through r and perpendicular to u. With this definition, P{f) 
is defined with respect to a uniform measure in the 2D space of possible forces. 
Given P{f,r), the local macroscopic stress tensor is given by PU] 

'Ta/3= / fdf d9no.n0fP{f,9,f). (2) 

Jo J-TT 

We wish to construct an equation that describes the ensemble average of the 
spatial variations in the chain densities for a system subject to specified boundary 
conditions. To do so, we consider the probability that a chain / will exist at a point 
f+ en, assuming we know P{f, fi, r). 

A given material will be characterized by a two scalar parameters, A and Y , and 
two angular functions, ips and 0/. 

• A is the average distance in any specified direction between the point where 
a chain begins and the nearest point that will cause it to split; i.e., the mean 
length chains would have if there were no fusions. 

• 4's ifi I /2, /a) is probability that a given sphtting of a /i-chain results in chains 
with strengths and directions given by /2 and /a, normalized to unity in the 
sense defined below. 

• 4'f{fi I f2, fs) is relative probability that a /2-chain and a /3-chain will fuse 
when they intersect and thereby form a /i-chain, normalized to unity in the 
sense defined below. 

• y is an overall efficiency with which two intersecting chains will fuse. That 
is, the total probability that two intersecting chains will fuse is Y(j)f 

(j)s and (j)f must both include delta functions that enforce force balance at each 
intersection and Heaviside functions that guarantee all forces are compressive. 
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Wc neglect, for now, the finite size of individual grains. That is, we assume that A 
is large compared to a grain diameter and treat the granular packing underlying the 
force chain network as a continuous medium. The equation governing the spatial 
variation of P as follows (for notational convenience, we have dropped the r from 
the argument of all of the P's): 



(n • V)P(/) = -J J fdf dO' f'df dO" \sm{0' - e")\ </.,(/ | /', f")P{f) (3) 

+\ j fdf d0' fdf de" \smie - e")\ Mf' I /, f")Pif') 
- 2Y j fdf de' fdf de" \smie - e") \ Mf I /, f")P{f)P{f") 
+ Y j fdf de' fdf de" |sin(e' - e") \ Mf\ f")P{f')P{f") 

Each term on the right hand side represents a type of event that can alter the 
density P{f). We discuss each in turn: 

• The first term represents the loss of /-chains due to splitting. The sin factor 
is exhibited explicitly for future convenience. Here P{f) may be taken out 
of the integral. To ensure that A is the mean distance between a chain will 
propagate before splitting (in the absence of interactions with other chains), 
we must therefore have 



/ 



fdf de' fdf de" |sin(0' - 0")\Uf\ f, /") = 1. (4) 



Note that (f)s has dimensions of Xj forced. 

The second term represents the creation of /-chains due to the splitting of 
chains in other directions. The splitting function here must be the same as in 
the first term, but with arguments exchanged, since the rate at which chains 
chains appear due to splitting must match the rate at which the parent chains 
split. The factor of 2 counts the identical integral arising from the exchange 
of the prime and double-prime labels. 

The third term, which is nonlinear in P, represents the loss of /-chains due to 
their fusions with /"-chains (or /'-chains, hence the factor of 2). The quantity 
|sin(^ — ^")| P(/)P(/") is the density of intersections of / and /"-chains, and 
(/)/ is specifies the probability of fusion when two such chains meet. Here again 
P(/) can be taken out of the integral. We choose to normalize (^j such that 
the remaining integral has a maximum value of unity over the set of P{f"ys 
of the form P{f") ~ 5{f" — f), for all values of f. In other words, we 
normalize 4>f according to the type of intersection most likely to produce a 
fusion if all chain densities were identical. In equation form: 

fdfd0'fdfd0"ct>jif\ff")S{f-fo) =1, (5) 

or, equivalently, 

^ = 1, (6) 



max 



Ifdfd0' Mf'lffo) 



where the maximization is over all possible values of / and /q. The param- 
eter Y then provides an absolute measure of the fusion efficiency. (j)f has 
dimensions 1/ force^, in contrast to (ps- 
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• The fourth term represents the creation of /-chains due to the fusion of chains 
in other directions. For consistency, the fusion function must be the same as 
in the third term. 

For materials composed of perfectly rigid grains, a rescaling of all of the forces 
in a given DFCN yields another perfectly acceptable network. Hence, Eq. © 
must be invariant under a rescaling of all the force intensities. This is verified by 
straightforward dimension counting (unlike the form of the Boltzmann equation 
suggested in Ref. 

Consider further the general forms of and (jjf. For an isotropic material, (j) 
can depend only on the differences between angles. The general form is 

0a(/ 1 /', /") - Faif, /', nsif + f" - muwmnMo' ~ &' ~ o), (?) 

where Fa is a combination of the force intensities that provides the correct dimen- 
sions for (j)a, and tpa is a symmetric function of its arguments. (One might imagine 
functions of the ratios of force intensities multiplying the arguments of "00, but 
these will always reduce to functions of the angles alone due to the S function.) 

In the case of fusions, dimension counting in Eq. implies Ff is dimensionless 
and hence equal to unity, up to a constant that can be absorbed in Y. In the case of 
splitting, on the other hand, dimension counting in Eq. I@J implies Fs has dimension 
1/ force^. The relations between the fs enforced by the S function guarantee that 
any combination of /, /', and /" with the right dimensions is equally valid; the 
differences between them can simply be absorbed into ipg. The natural choice is 
Fs = !/(/'/")• With this choice a constant i/^s — I/tt^ corresponds to an equal 
probability for every possible splitting configuration. 

It is clear that A can be scaled to unity without loss of generality by choice of the 
unit of length. Though 1^ is a dimensionless parameter with physical significance, 
from a mathematical point of view it plays a trivial role in Eq. (PJ . To solve for P 
for any given value of Y, we solve the case Y = 1, then simply divide all P's by Y. 
From here on, we take A = 1 and Y = 1. 

2.3. Specialization to discrete directions. We now make two simplifying as- 
sumptions to arrive at a set of ordinary differential equations that permits analytical 
solution. First, we assume that and 0/ are nonzero only for vertices of the forms 
shown in Figure [31 This means that all chains lie in the discrete set of directions 
9n = {n — ^)'K I {2N) (measured from the positive z direction, which is downward), 
and that the strong force at a given vertex is always related to two weak ones by 
a factor of ^ = 2cos(a). Defining = fo£,"\ the continuous function P{f,9) be- 
comes a set of discrete functions P{fm-, 6n), which will be denoted Pn,mSif — fmfi-n)- 
In this notation, it is always assumed that n is taken modulo AN. Note that Pn,m 
is simply a number per unit length; the dimensions of force are taken care of by the 
two-dimensional delta function. We will treat the case a = 7t/{2N) explicitly here. 
Generalization to any a that is an integer multiple of this is straightforward. 

This first assumption corresponds to the following forms for the angular parts 
of the splitting and fusion functions: 

MV'.V") = ljS{<p'-a)S{v" + a) + Hv' + a)Hv"-a)]. (8) 
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Figure 3. Star of chain directions and allowed vertices for solvable 
models. The star has AN-iold symmetry. The vertices show a 
splitting and a fusion with angle a — 7t/{2N). 

The difference in character between 4>s and 4> / stems from the fact that 4> f appears 
in integrals that are quadratic in P. When P is a sum of two-dimensional delta func- 
tions, the fusion integrals have two more delta function factors than the splitting 
integrals. The same can be seen in the normalization conditions of Equations 
and Equations ©• 

Second, we assume that the boundary conditions of interest are uniform across 
the top and bottom of the slab. Translational symmetry in the x-direction then 
dictates that all horizontal gradients of Pn^m vanish. 

Inserting these assumptions into Eq. Q and integrating both sides over a small 
volume of force space in the vicinity of fmnn, we obtain the following ordinary 
differential equations for the Pn.m's: 

^Pn,m + -P)i+l,?Ti+l + Pn-l,m+l (10) 

The result is straightforward, but one must be careful to account for all the trigono- 
metric factors arising from integrating over the delta functions, including the inte- 
gral on both sides that is required to isolate Pn,m- 

The different terms on the right hand side of Eq. ifTIHl correspond to the processes 
described above that can create or destroy chains. In the present case, the first term 
represents the decay of Pn,m along the direction 0n due to the splitting of an existing 
(n, m) chain. The next two (linear) terms account for the addition to Pn,m due to 
the splitting of other chains in adjacent orientations. Note that if a chain is to split 
and create a contribution to Pn,m, which must have strength /oC™j it must begin 
with strength /o'C™"''^- The first nonlinear term accounts for events in which two 
chains of strength /oC™~^ fuse, and the last two nonlinear terms to events in which 
chains contributing to Pn.m are deflected due to fusions with other chains. 

In Reference Eq. was studied in detail for the case of 6- fold symmetry. 
In that case we have ^ = 1, so the hierarchy of equations indexed by m collapses. 
In addition, cos 02 = cos 6*5 — 0, so the entire system reduces to four ODEs and two 
algebraic relations. If we assume mirror symmetry about the vertical axis, these 
are reduced to two equation that can be solved completely. 

In order to observe the variations of P with force strength, we must consider 
symmetries other than 6-fold. To avoid technical complications arising from the 
presence of strictly horizontal chains having cos 6 = 0, we consider only the cases 
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of 4iV-fold symmetry for positive integers N. The simplest case is that of 8-fold 
symmetry, and example of which is illustrated in Figure ^ 

Note that the discreteness of the possible orientations of force chains docs not 
imply any discreteness in the possible positions of the splitting or fusion vertices. 
Note also that the system remains isotropic in the sense that Eq. ^| including all 
coefficients, is identical for all n, except for the geometric factor on the left-hand 
side that results from restricting attention to variations in z and not x. 

3. Homogeneous networks for symmetric splittings and fusions. General 
solutions of Eq. (|10|) are not yet known. We study here the important special 
case of homogeneous solutions; i.e., the fixed points, for which all gradients vanish. 
These solutions are both nontrivial and crucial for understanding the response 
function. Consider the trajectory of the vector Pn,m as z is varied. As detailed in 
Reference 21 for the 6-fold case, trajectories that do not pass very close to a stable 
(or marginally stable) fixed point inevitably lead to divergences at z >> A that are 
reflected in unphysical negative values of some Pn,m- Trajectories that avoid these 
divergences do so by staying in the vicinity of a fixed point over most of the range 
of z. As also noted in Reference [21], the trivial fixed point Pn,m = is unstable 
and hence cannot serve as a reference state for a linear response theory. Thus the 
response function we ultimately seek will be dominated by the DFCN structure 
corresponding to some nontrivial fixed point. 

3.1. Isotropic solutions. It is instructive to begin with a search for isotropic 
solutions; i.e., solutions for which P„_m — Pm for all n. The fixed point equation 
derived from Eq. (|10|l reads 

-P„ + 2P„,+1+P2_i-2P2 =0. (11) 

To simplify this recursion relation, define 

r/ — P — 

Eq. pil) says 2dm+i — dm = 0, which implies 

dm = dt:,2 ™. 

For do ^ 0, Eq. lO gives 

Pm+l = Pm + da"^ ^ 

For do < 0, this recursion relation either produces a negative numbers or a diver- 
gence at large m, neither of which is physically acceptable. If Pq is sufficiently 
small, the repeated squaring drives Pm toward zero at least as fast as with 
p > 1. The negative contribution from the do term can only speed up this decay 
as long as Pm is positive. But this means that Pm decays faster than 2~™, so the 
do term eventually makes P„i negative. Now if Pg is large enough to avoid Pm 
decaying toward 0, repetitive squaring causes it to diverge like with p > 1; the 
do term becomes irrelevant at large m. 

For do > 0, on the other hand, problems arise for large, negative m. We write 
Pm-i = VPm ~ do2~'" and see that complex values will be generated unless Pm 
is always greater than do2~'". But this is impossible because Pm-i is strictly less 
than Pm- 

Thus we are left with do = as the only physically relevant case. When do = 0, 
Eq. (|l()|l is satisfied in a special way: both —Pm + Pm-i = and 2P,„+i — 2P,^ = 



(12) 
(13) 
(14) 
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are satisfied simultaneously. These two relations are actually identical up to a shift 
in TO, and admit a one-parameter family of solutions 

P^=p-^'\ (15) 

where we must have p > 1 to avoid divergence at large to. Despite its apparent 
simplicity, it is instructive to examine this solution in a bit more detail. We consider 
first the behavior at large negative m, then the behavior at large positive to. 

Note that Pm 1 for large negative to. This means that the densities of all 
chains of intensity fo^~^ for large to are the same. In other words, the total 
density of force chains, Pm diverges. Though this may appear troublesome, 
it does not yield a divergence in the pressure: Pmfm remains perfectly finite. 
The divergence in chain density is an artifact of the restriction to a fusion function 
(pf that does not allow weak chains to fuse with stronger ones. Removal of this 
artifact requires working with a continuous distribution of force chain directions 
and intensities and it beyond the scope of this work. 

For large positive m, we can compute the probability distribution for contact 
forces, a quantity that has been shown in a variety of experiments and numerical 
simulations to decay exponentially for large forces. For the homogeneous, isotropic 
solution, we have 

V{fm)^ Pm/ifm- fm^l), (16) 

where fm = /oC™ and Pm = ■ The denominator in this expression is just the 
spacing between points with successive m's, the required conversion factor from the 
discrete density Pm to the density per unit force. This can be rewritten as 

V{f)^p-U/fof/f^ with (17) 

Recall that p > 1 and ^ < 2, so /3 > 1 and V decays faster than exponentially for 
large /. For the 8-fold case, we have (3 = 2 and hence a Gaussian decay. But smaller 
splitting angles give ^ closer to 2, and for a splitting angle as large as 2q; = 60°, we 
have (3 ~ 1.26. In this case, the full distribution (including the 1// term) can look 
surprisingly like a simple exponential over several decades, as shown in Figure 0] 

Experience with the 6-fold model suggests that the isotropic fixed points of 
Eq. (|10|l are non-generic. In the 6-fold model, linearization of the theory in the 
vicinity of a fixed point generally leads to hyperbolic response functions with peaks 
that are sharper for more strongly anisotropic fixed points. The isotropic case is 
pathological in that the peak width is divergent [21] and additional cancellations 
conspire to produce single peak, but one that obeys different scaling laws from the 
predictions of elliptic models. [21] It is therefore important to consider solutions of 
Eq. H10|) that are not constrained to be isotropic. 

3.2. Anisotropic solutions. To find fixed points of Eq. I|10|l . we begin with the 
following ansatz : 

P„,™=P«'"+""'^", (18) 

where p is a real parameter and qm, Xn, and <^ are to be determined. (This guess 
combines features of the isotropic solution above and the 6-fold solutions discussed 
in Reference [H].) Substituting into Eq. (|10|l . we have 
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Figure 4. Single contact force distributions for 8-fold and 12-fold 
models. For each case, curves are shown for p — 1.2, 2, and 5. The 
straight, dashed lines are guides to the eye to emphasize the fact 
that the curves look very close to ordinary exponential decays. 



Taking a cue from the isotropic case, we consider the possibility that the cancellation 
required by this equation occurs term by term in the following way:'^ 

p9m+2:„C™ _ p2q„_i + (rE„+i+rE„_i)C""\ J^20) 
pgm+l+a:„_iC'"+^ _ p2qm + {x„+x„-2)C" (^21) 
pq,„+i+a;„ + iC'"+^ _ p2q„^ + {x„+x„+2)C"^ ^ ("22) 

Note that all three of these equations are identical up to shifts of n and m. Since 
the q terms are independent of n, we obtain two recursion relations: 

g„i+i = 2g„; (23) 

XnC = Xn-l+Xn+l- (24) 

(25) 

Eq. 1221) implies 

qra = -qo^^, (26) 
where the negative sign defines a convention for the sign of go- Eq. 1)24(1 is an 
eigenvalue equation for a AN x 4iV matrix M . The eigenvalues Q and eigenvectors 
yi"'^ of M are easily found: 

4. -^cos(^2^J , j = 2N+l,..AN-l ^^'> 

Thus for each of the AN eigenvalues and eigenvectors, wc have what appears to be 
a three-parameter family of fixed point solutions of Eq. (|10() of the form 

P„,™=p-*''"+"°""'''^^'", (28) 
the free parameters being p, go, and xq. Recall that the force intensity associated 
with the chains whose density is Pn,m is fm — fo€^ and note that ^ = Ci- Note 



•^Unlike in the isotropic case, we do not prove here that all physically plausible solutions must 
have this form. 
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also that for j — 0, we have C = 2, so that the xq and go terms can be combined. 
This leaves us with precisely the same isotropic solution discussed above. 

Eq. H28(l requires further examination to determine the range of different struc- 
tures it represents. First, each three-parameter family is really only a one-parameter 
family of distinct physical networks. There are ways of rescaling the parameters 
that have no effect on the physical system begin described: 



The first of these corresponds to a redefinition of the force scale, /o — > /oC'^, which 
cannot affect the physics. Strictly speaking, the system is exactly invariant under 
this transformation only if is an integer. Arbitrary values of v, however, merely 
shift the positions of the discrete set of P„,m's along a single smooth curve for each 
n, with integral i^'s being the values that shift the entire set into itself. For our pur- 
poses, the differences between solutions related by non-integer are unimportant 
details. (See FigureEl) The second transformation leads to a mathematically iden- 
tical solution. (The case go = is not physically relevant, as it leads to divergence 
of the total stress from large m force chains.) By performing the first transforma- 
tion with V = log(a;o/<Zo)/ log(2/C), followed by the second, we can adjust go and 
Xq both to unity. 

Second, note that we cannot take arbitrary linear combinations of the solutions 
to Eq. (|24|l . since they would not satisfy the nonlinear Eq. HlOf) . We can, however, 
take linear combinations of eigenvectors that have degenerate eigenvalues. Since 
CiN-j = Cj foi" 1 ^ J ^ ~ 1 can generate solutions from arbitrary linear 
combinations of the degenerate pair of eigenvectors y^^^ and y^'^^^iK Thus j = 
and j = 2N each provide a single continuous family of solutions parameterized by 
p, while all other j's come in pairs that provide families parameterized by p and by 
the relative weights of the two if-^^^s. 

Consider now the family spanned by y^^' and for some specific j. By 

forming linear combinations, we can construct basis vectors ej that are symmetric or 
antisymmetric under reflection through a vertical axis. The symmetric combination 
will give solutions for which Pn,m = P'iN+i-n,rm as might be expected if the average 
confining force on the plates is purely normal. The antisymmetric combination will 
give solutions that break this symmetry, though they do not yield any negative 
values of Pn,m- These might be important for describing systems subjected to 
shear as well as compression. The explicit forms of the symmetric combinations 
(including the trivially symmetric isotropic case) can be expressed as 



a stress configuration with a dipole-like anisotropy. At each generation m, the 
angular variation Pn,m has a single maximum (for p > \) aX n — 1 and a minimum 
n = 2N. There are more downward-pointing than upward-pointing chains for 



(90, 2^0 ) 

{p,qo,xo) 



(2''go,r.To) and 
(p«,l,a;o/go). 



(29) 
(30) 




1. 




Figure 5. Anisotropic solutions for 8- fold case. The solutions pic- 
tured have reflection symmetry about the vertical; i.e., Pn.9-m — 
Pn,m- Left: Two j = 1 "dipole-like" solutions. Black circles show 
a solution corresponding to p = 10, Xq = qo = 1. Grey circles show 
a solution for xq — 1.7, with p adjusted and m shifted to show 
that it belongs in the same family as the = 1 solution. Right: 
A solution with a higher order anisotropy, j — 3. The oscillation 
of Pn,m with TO is apparent. The solid line is a guide to the eye for 
following the oscillations of Pi^m- 



each undirected orientation. Moreover, because Ci is positive, the decay of Pn^m 
with TO is monotonic both at large to, where it decays to 0, and at large negative 
TO, where it decays to 1. Figure [S] shows the form of Pn,m for this case. Note 
that the term "dipole-like" applies here to the P's. Since the stress field does not 
distinguish between force chains in opposite directions, the anisotropy in the stress 
is "quadrupole-like" - stronger in the vertical directions than in the horizontal. 

A striking feature of the solutions is the occurrence of peaks for some n, which 
become more pronounced for larger p. Peaks are evident for all n when the stress, 
fmPn.m, is plotted rather than the chain density. (See Figure 0) Differentiating 
Eq. 1)28(1 with qo = xq = 1, we find that the peak in Pn,m occurs at 

fogyn ^ + loglog |C„| - loglog2 

log2-log|C„| ^ ^ 

which is independent of p. (The peak in fmPn,m does depend weakly on p.) The 
force scale represented by this peak must have its origin in the boundary conditions, 
since a linear rescaling of all forces in the bulk has no effect on the DFCN structure. 
Recall that part of the process required to collapse the solutions onto the one- 
parameter family indexed by p was to rescale /o, which corresponds to a rescaling 
of the strengths of the force chains injected at the boundaries. There are subtleties 
lurking here, however, as will be discussed further in Section 

A full calculation of the single contact force distribution V{f) is slightly more 
complicated than for the isotropic case, as it requires a sum over the different chain 
orientations. It is not hard to extract the dominant term at large /, however. For 
the symmetric j — 1 case with qo — xq ~ 1, Eq. (|28|) can be written as 

P„,„=p-/t+^i^>/^, (34) 

where (3 = log 2/ log ^ as above, and 7 = log (^1/ log ^ = 1. Since f3 > 1, the 
dominant term at large / will be Pi^m — p^f^^^f . 




Figure 6. Stresses for an anisotropic solution. The solution pic- 
tured is the same as the one on the left in Figure El Here we 
show the compressive stress associated with a given chain direc- 
tion rather than just the chain density, plotted vs. m and fm. The 
dotted lines on the right are just guides to the eye. 



It is important to note that the anisotropy of the solution is not a consequence of 
any material anisotropy. The bulk properties of the material we are describing are 
perfectly isotropic. The boundary conditions on the system will generally not be 
isotropic, however, so there is no reason to expect that the pre-stressed state about 
which the response function should ultimately be calculated should have isotropic 
stress. The generic situation described by the discrete DFCN theory presented 
here is one in which an intrinsically isotropic material is subjected to loading forces 
that generate anisotropic distributions of stress chains. Because different degrees of 
anisotropy correspond to different fixed points of the full nonlinear theory and hence 
to different linearizations, the pre-stressed state may have a profound influence on 
the response function. This is in contrast to an ordinary elastic medium, in which 
the response function is not affected by an anisotropic background stress. 

The solutions for higher j have features that seem a bit unlikely from a physical 
standpoint, though not obviously unacceptable. For 2 < j < iV, we have higher 
orders of anisotropy, with P„^rn exhibiting additional maxima and minima as a 
function of n at each m, but still with monotonic decay at large to. For N < j < 2N, 
the eigenvalues are negative, which causes C"* and hence Pn,m to oscillate as to 
varies. Whether or not these solutions are relevant hinges on the question of which 
boundary conditions are associated with generic physical systems. 

3.3. Remark on boundary conditions. Given several families of homogeneous 
solutions, it is natural to ask which solution will be selected for a given set of 
boundary conditions. From a mathematical perspective, this is a straightforward 
question. As described in detail for the 6-fold case in Reference , the only way to 
specify boundary conditions on Eq. (fTn|l that are sure to be self-consistent and not 
produce negative chain densities is to specify the densities of only the ingoing chains 
at each boundary. In general, the densities specified will not match perfectly any 
of the homogeneous solutions Pn,m for all m and n, even if horizontal translational 
symmetry is assumed. Based on experience with the 6-fold case, we might expect 
the trajectory in the space of Pn,m as the slab is traversed to quickly approach a 
fixed point and stay near it over a substantial portion of the slab, then make a 
rapid transition to the vicinity of another fixed point, and finally move off to a 
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point consistent with the bottom boundary eondition. The details of the transition 
regions and selection of the fixed point remain to be worked out in the present case. 

From a physical perspective, the situation is not so clear. Consider the following 
two possible choices of boundary conditions for the top surface. In both cases, 
assume that Pn,m = for all ingoing {n,m) except the ones specified as follows: 
(BCl) Po,5 = -Pi,5 = 1; and (BC2) Pq.i = -Pi.i = C^- Assume also that the boundary 
conditions on the bottom slab are simple reflections of these through the horizontal. 
Thus we have boundary conditions with the same total stress stipulated, but in one 
case it is injected via a sparse set of strong forces, and in the other via a denser set 
of weaker forces. These conditions will not necessarily lead to the same fixed points 
in the bulk. Since the chain densities are measured with respect to the intrinsic 
length scale A, BCl and BC2 are not simply related by a scale transformation. We 
therefore expect them to lead to fixed points corresponding to different p's, which 
implies that they will require different rescalings of /q. This in turn indicates that 
the force scales fixed by the two different boundary conditions differ, in spite of the 
fact that they correspond to specifying the same overall force. 

The resolution of this apparent paradox lies in the fact that the boundary condi- 
tions do not specify the densities of the outgoing chains. Ultimately, the force scale 
must be determined by the zz component of the stress, which must be constant 
throughout the slab; i.e., must determine the prescribed rescaling of /q. But 
the densities of outgoing chains produced by BCl and BC2 will differ, which will 
cause (Jzz to differ in the two cases. In order to determine cr^z, which is the total 
pressure applied to the top surface, we must first determine the full solution, or at 
least which fixed point the trajectory approaches. 

The converse situation, in which the stress is specified but not the individual 
chain densities, is more familiar, but just as subtle. One might expect it to be 
sufficient to specify some components of the stress tensor and its spatial derivatives 
on the top boundary and others on the bottom, but the DFCN equations require 
independent specifications of densities of chains of many directions and strengths. 
In an experiment in which a granular material is put under compression between 
to fiat plates, it is not at all clear how the chains densities at the boundaries are 
determined. Both the distribution over directions n and strengths m are relevant. 
An important topic of future research will be to develop a method of assigning 
directions to force chains observed in experimental images or numerical models, 
and to learn from them how the boundary conditions appear to be imposed. 



4. Conclusion. In this paper we have exhibited for the first time solutions of the 
Boltzmann equation for the chain densities in a model that permits the interaction 
of chains of different intensities. These solutions have highly nontrivial structure 
that raises a number of interesting questions. Most importantly, they permit cal- 
culations of both stresses and single particle force distributions in a unified way. 
Refinement of the model is clearly necessary: wc would like to avoid divergences in 
the density of weak chains; and we would like to generalize the splitting and fusion 
functions, which might bring the single contact force distribution V{f) into closer 
agreement with experiment. Nevertheless, there is much to be gained from further 
study of the models presented here, or minor modifications of them. 

The availability of closed form analytic solutions suggests that a complete the- 
ory of the stress configurations and response functions for these particular models 
is within reach. In addition, it appears possible that the linearized theory in the 
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vicinity of the fixed points presented here could lead to a new derivation of consti- 
tutive relations for closing the equilibrium equations for the macroscopic stresses. 
Finally, the solutions presented here may provide a basis for extensions to similar 
discrete models in which more than one type of vertex structure is allowed. All of 
these topics are currently under study. [53] The author looks forward to continued 
collaboration with David Schaeffer on these topics, as he is certain to offer valuable 
insights and pose questions that lead to fruitful calculations. 
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